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Abstract. We study long-range interacting systems perturbed by external stochastic 
forces. Unlike the case of short-range systems, where stochastic forces usually act locally 
on each particle, here we consider perturbations by external stochastic fields. The system 
' ^ I reaches stationary states where external forces balance dissipation on average. These states 

do not respect detailed balance and support non-vanishing fluxes of conserved quantities. 
. We generalize the kinetic theory of isolated long-range systems to describe the dynamics 

of this non-equilibrium problem. The kinetic equation that we obtain applies to plasmas, 

■ self-gravitating systems, and to a broad class of other systems. Our theoretical results hold 
^ ! for homogeneous states, but may also be generalized to apply to inhomogeneous states. We 

obtain an excellent agreement between our theoretical predictions and numerical simulations. 

■ We discuss possible applications to describe non-equilibrium phase transitions. 

PACS numbers: 05.20.Dd, 05.70.Ln, 05.40. 
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1. Introduction 

Most physical systems are out of equilibrium either because of couphng to thermal 
baths at different temperatures or because of external forces that break detailed balance. 
Studying non-equilibrium stationary states is an active area of research in modern statistical 
mechanics. It is indeed a lasting challenge to achieve for non-equilibrium systems a level of 
theoretical understanding similar to the one established for equilibrium systems [THS]- 

In this Letter, we consider systems of particles interacting through two-body non- 
integrable potentials, also called long-range interactions. Examples include plasmas and self- 
gravitating systems (globular clusters, galaxies), where particles interact through repulsive 
or attractive Coulomb and attractive Newton potential, respectively. Our work also applies 
to a large class of models with non-integrable interactions, such as spins, vortices in two 
dimensions, and others, which have been studied extensively in recent years [IHH]. 

Systems with non-integrable potentials are often forced through external stochastic 
fields. For example, globular clusters are affected by the gravitational potential of their 
galaxy, thereby producing a force that fluctuates along their physical trajectories. In 
addition, galaxies themselves feel the random potential of other surrounding galaxies, and 
their halos are subjected to transient and periodic perturbations due, for example, to the 
passing of dwarfs or to orbital decaying [9]. Plasmas may also be subjected to fluctuating 
interactions imposed by environmental electric or magnetic fields [10]. These physical 
situations often lead to a stationary state where the power injected by the external random 
fields balances on average the dissipation. To the best of our knowledge, such non-equilibrium 
stationary states in systems with non-integrable potentials have not been studied before, and 
this work provides a first step in this direction. 

Unlike systems with short-range interactions, stochastic perturbations in long-range 
interacting systems often act coherently on all particles and not independently on each 
particle. Moreover, unlike short-range systems, it is not natural to consider long-range 
systems as being coupled to thermal baths at the boundaries. Thus, the non-equilibrium 
stationary states that we study are rather different from the ones in systems with short- 
range interactions. These states do not verify detailed balance and support non-zero fluxes 
of conserved quantities, which are basic ingredients of non-equilibrium stationary states. 

Theoretical results on isolated systems with long-range interactions include the kinetic 
theory description of relaxation towards equilibrium. In plasma physics, this approach 
leads to the Lenard-Balescu equation or to the approximate Landau equation [HI [12]. 
These equations, or some of their approximations, are grouped as the collisional Boltzmann 
equation in the astrophysical context. The main theoretical result of this Letter is a 
generalization of the kinetic theory to describe non-equilibrium stationary states, valid for 
small external perturbations and spatially homogeneous stationary states. 

The non-equilibrium kinetic equation that we obtain describes the temporal evolution 
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of the one-particle distribution function. When the system is not far from equihbrium, it is 
natural to expect that the system settles into a stationary state. We find that in such a state, 
the one-particle momentum distribution is non-Gaussian. The kinetic equation describes the 
evolution of the kinetic energy, and its prediction of the stationary state compares very well 
with A^-body numerical simulations. 



2. Stochastically forced long-range interacting systems 

Consider a system of particles interacting through a long-ranged pair potential. The 
Hamiltonian of the system is 

N 2 1 ^ 

where qi and pi are, respectively, the coordinate and the momentum of the z-th particle, 
while v{q) is the two-body interaction potential. The particles are taken to be of unit mass. 
In this paper, for simplicity, we consider g^'s to be scalar periodic variables of period 27r; 
generalization to qi G M", with n = 1, 2 or 3, is straightforward. 

In self-gravitating systems, since the dynamics is dominated by collective effects, it is 
natural and usual to rescale time in such a way that the parameter 1/N multiplies the 
interaction potential [13]. In plasma physics, the typical number of particles with which one 
particle interacts is given by the coupling parameter F = nA|), where n is the number density 
and Xr, is the Debye length. It is then usual to rescale the time such that the inverse of a 
power of r multiplies the interaction term [H]. These reasons justify the rescaling of the 
potential energy by 1/N in Eq. ([1]), known as the Kac scaling in systems with long-range 
interactions |14j . 

We perturb the system ([T]) by the action of the stochastic field F{qi,t). The resulting 

equations of motion are 

OH ^ OH X 

gi = — , and Pi = api + y/a F{qi,t), (2) 

Opi oqi 

where a is the friction constant, and F{q,t) is a statistically homogeneous Gaussian process 

with zero mean and variance given by 

{F{q,t)Fiq',t')) = C{\q-q'm-t'). (3) 

The hypothesis that the Gaussian fields are statistically homogeneous, i.e., the correlation 
function depends solely on \q — q'\, holds for any perturbation which does not break space 
homogeneity. Such a hypothesis will also be essential for the following discussions where we 
consider homogeneous stationary states. Now, C{q) represents correlation, and is therefore 
a positive-definite function |T5j. Its Fourier components are thus positive: 

Ck = — dq Ciq)e~'''' > 0, = Cq + 2 V cos(A;g). (4) 



It will be convenient to use the following equivalent Fourier representation of the Gaussian 
field F{q,t): 

oo 

F{q, t) = y^Xo + J2 [cos(A;g) + sm{kq) Y^] , (5) 

k=l 

where Xk and Yk are independent scalar Gaussian white noises satisfying {Xk{t) Xk'{t')) = 
5k,k'5{t - t'), {Yk{t) Ykit')) = 5k,k'5{t - t'), and (Xfc(t) = 0. 

Using the Ito formula [16] to compute the time derivative of the energy density e = H/N 
and averaging over noise realizations give 

where k = YliLiPi /C^^) kinetic energy density. The average kinetic energy density 

in the stationary state is thus (k)^^ = C(0)/4. 

In the dynamics ([2]), fluctuations of intensive observables due to stochastic forcing are 
of order -y/a, while those due to finite-size effects are of order 1 / \/N. Moreover, the typical 
timescale associated with the effect of stochastic forces is 1/a (see, e.g., Eq. ([6])), while the 
one associated with relaxation to equilibrium due to finite-size effects is of order A^, see [ll[5] . 

In the following, we analyze the dynamics ([2]) in the joint limit — t- oo and a — 0. 
While the first limit is physically motivated on grounds that most long-range systems indeed 
contain a large number of particles, the second one allows us to study non-equilibrium 
stationary states for small external forcing. Moreover, for small a, we will be able to develop 
a complete kinetic theory for the dynamics. 

For simplicity, we discuss in this Letter the continuum limit Na ^ 1, when stochastic 
effects are predominant with respect to finite-size effects. Generalization to other cases {Na 
of order one, or, Na <^ 1) is straightforward, as discussed in the conclusion. 



3. Kinetic theory 



A natural framework to study the dynamics ([2]) is the kinetic theory. We now describe 
the theoretical approach to derive this theory, while some of the technical results will be 
explicitly obtained in a longer paper [17]. The central result is the kinetic equation ( fTOl) 
below, which describes the time evolution of the single-particle distribution function. 

We consider the Fokker-Planck equation associated with the equations of motion ([2]). 
It describes the evolution of the A^-particle distribution function fN{qi,---,qN,Pi,---,PN,t) 
(after averaging over the noise realization, Jn is the probability density to observe the system 
with coordinates and momenta around the values {qi.,Pi}i<i<N at time t). This equation can 
be derived by standard methods [16]. We get 

N r 
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We have proved by analyzing the so-called potential conditions [18] for this Fokker-Planck 
equation that a sufficient condition for the stochastic process ([2]) to verify detailed balance 
is that the Gaussian noise is white in space, that is, Ck = c for all k. This condition is not 
satisfied for a generic correlation function C. Steady states are then true non-equilibrium 
ones, with non-vanishing currents and a balance between external forces and dissipation. 

Similar to the Liouville equation for Hamiltonian systems, the A^-particle Fokker- 
Planck equation is a very detailed description of the system. Using kinetic theory, 
we want to describe the evolution of the one-particle distribution function f{zi,t) = 
lY[iL2 fN{zi,--.,ZN,t) (we use the notation Zi = {qi,Pi) whenever convenient). Note 
that the normalization is J dz f{z,t) = 1. 

In plasmas and self-gravitating systems, due to the long-range nature of the interactions, 
the one-particle distribution function is not affected by the two-particle distribution function 
at leading order in and therefore, its evolution is described at leading order by the 

Vlasov equation. Finite-size effects however induce weak correlations whose effects on the 
long-time evolution of the one-particle distribution can be computed self-consistently in the 
framework of the kinetic theory by using perturbation theory. A complete treatment of the 
problem leads to the Lenard-Balescu equation [II1II2]- In a similar way, for our problem, the 
evolution will be described at leading order by the Vlasov equations due to the long-range 
nature of the interactions. Weak stochastic forces lead to weak correlations that affect the 
long-time evolution. This case can be treated by following a generalized kinetic approach, 
as we now describe. 

Substituting in the A^-particle Fokker-Planck equation (JTj) the reduced distribution 
function fs{zi,...,Zs,t) = J YliLs+i fN^Zi, zj^,t), and using standard techniques [19], 
we get a hierarchy of equations, similar to the BBGKY hierarchy. We split the 
reduced distribution functions into connected and non-connected parts, e.g., f2{zi,Z2,t) = 
f{zi, t)f{z2, t) + ag{zi, Z2,t), and then neglect the effect of the connected part of the three- 
particle correlation on the evolution of the two-particle correlation function. This scheme 
is consistent at leading order in the small parameter a, and is the simplest closure scheme 
for the hierarchy. For simplicity, we moreover assume that the system is homogeneous: / 
depends on p, and g depends on \qi — q2\, Pi and p2, only. The ffist two equations of the 
hierarchy are then 



df d a d^f d f 

a— J dqdp2 v'{q)g{q,p,p2,t), 

and 



dt dp 2 dp' 

j dgsdps v'{qi - q^ijgiqz - q2,P3,P2,t) 



dg_ 
dt 



dg df 
dqi dp 
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Ciqi - q2] 



dl 
dp 



pi 



dp 



(9) 
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where the symbol {1 ^ 2} means an expression obtained from the bracketed one by 
exchanging 1 and 2, while the prime denotes differentiation. 

To obtain from these equations a single kinetic equation for the distribution function 
/, we have to solve Eq. (|9]) for as a function of / and plug the result into the right hand 
side of Eq. From these two equations, we readily see that the two-particle correlation g 
evolves over a timescale of order one, whereas the one-particle distribution function fXp,t) 
evolves over a timescale of order 1/a. We use this timescale separation, and compute the 
long-time limit of g from Eq. iQ by assuming / to be constant. This procedure is equivalent 
to making the Bogoliubov's hypothesis for deriving the kinetic theory of isolated long-range 
systems. For the timescale separation to be valid, it is also required that the one-particle 
distribution function /(p, t) is a stable solution of the Vlasov equation at all times. 

The solution of equations of the type ([9]) is quite technical (see the long appendix in 
Nicholson's book [11]). Equation iQ differs from the corresponding equation for an isolated 
long-range system in that the term on the right hand side is different in the two cases, and 
cannot be solved by methods known in the literature. The main technical achievement that 
aided this work is to be able to solve Eq. ([9]). In a future paper [17], we will give the details 
on the solving procedure. In brief, the method relies on making a parallel between the 
Lyapunov equations for infinite- dimensional Ornstein-Uhlenbeck processes and their general 
solutions, and Eq. (Q. Using this method, we get the desired kinetic equation: 
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Here, is the k-th Fourier coefficient of the pair potential v{q), the quantity is defined 
in Eq. (jl]), while J* indicates the Cauchy integral, and the dielectric function e is 

1 df 



e{k, u) 



lim 



1 — 2iiivkk i dp 



—i{uj + 17]) + ikp dp 



(12) 



The kinetic equation f fTOj) is the central result of the Letter. 

This kinetic equation has the form of a non-linear Fokker-Planck equation, since the 
diffusion coefficient D[f]{p) itself is a function of the unknown distribution function /. As 
Eq. (fTT]) shows, this coefficient has two parts, namely, (i) a linear part, C(0)/2, which is due 
to the mean-field effect of the stochastic forces, and (ii) a non-linear part due to correlations 
induced in the system by the stochastic forces. The contributions of different modes of the 
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stochastic force are independent of each other, that is, contributions proportional to do 
not couple with Vk' with k ^ k' . 

For consistency, the prediction of the evolution of the kinetic energy from the kinetic 
equation has to agree with Eq. ([6]). We have checked this by using Eqs. ( ITOj) and ( ITT]) , 
and proving that the integrals in the non-linear part of the diffusion coefficient give no 
contribution to the kinetic energy. 

As already mentioned, and is evident from Eq. ffTOj) . the time scale for the kinetic 
evolution is 1/a. This has been checked by performing direct numerical simulations, see 
Fig. [T] and the next section. Thus, a can be eliminated from the kinetic equation f lTU]) by 
a redefinition of time. Therefore, even for vanishingly small value of a, if a stationary 
distribution exists, it will be at distance of order one from a Gaussian momentum 
distribution. 




0.1 1 10 0.1 1 10 

at at 



Figure 1. (a) Kinetic energy density (k) and (b) {p^) as a function of at, for the values 
C(0) = 1.5 and ci =0.75. The data for different N and a values are obtained from numerical 
simulations of the stochastically forced HMF model, and involve averaging over 50 histories 
for A'' = 10^ and 10"^ histories for N — 10^. The data collapse implies that a is the timcscalc 
of relaxation to the stationary state. The inset shows the data without time rescaling by a. 

While a linear Fokker-Planck equation with non-degenerate diffusion coefficient can be 
proven to converge to a unique stationary distribution [18], this is not true in general for 
non- linear Fokker-Planck equations such as Eq. (fTOj) . We expect that if the system is not 
too far from equilibrium, the kinetic equation will have a unique stationary state. Far from 
equilibrium, the kinetic equation could lead to very interesting dynamical phenomena, such 
as bistability, limit cycle or more complex behaviors. The main issue is then the analysis of 
the evolution of the kinetic equation. Although some methods to study this type of equations 
exist [20], in order to provide some preliminary answers, we have devised a numerical iterative 
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scheme to compute some of the stationary states of the kinetic equation (fTOj) . We now 
describe the scheme. 

A hnear Fokker-Planck equation whose diffusion coefficient D{p) is strictly positive 
admits a unique stationary state 



For a given distribution fn{p), we compute the diffusion coefficient Dn{p) from Eq. (fTTl) . and 
then fn+i using Dn and Eq. (fT^ . This procedure defines an iterative scheme. Whenever 
convergent, this scheme leads to a stationary state of Eq. ffTOj) . Each iteration involves 
integrations, so we expect the method to be robust enough when starting not too far from 
an actual stationary state. However, we have no detailed mathematical analysis yet. 

In the next section, we discuss numerical results on A^-particle simulations, and the 
computation of stationary states from the iterative method mentioned above. 

4. Stochastically forced HMF model 

Until now, we have presented our theoretical analysis for a general two-particle interaction 
v{q). In order to perform simple numerical simulations, we now consider the case of the 
stochastically forced attractive Hamiltonian mean-field (HMF) model, which corresponds to 
the choice v{q) = 1 — cosg. 

The HMF model serves as a paradigm to study long-range interacting systems, and 
describes particles moving on a circle under deterministic Hamiltonian dynamics [2T||22]. 
This model has been studied a lot in recent times. It displays many features of generic 
long-range interacting systems, such as the existence of quasistationary states [22]. 
In equilibrium, the system displays a second-order phase transition from a high-energy 
homogeneous phase to a low-energy inhomogeneous phase at the energy density Cc = 3/4. 

Since the Fourier transform of the HMF interparticle potential is, for k ^ 0, Vk = 
— [Sk,i + /2, where 6k,i is the Kronecker delta, we see from the kinetic equation (ITOll that 
only the stochastic force with wave number k = 1 contributes to the non-linear part of the 
diffusion coefficient; all the other stochastic forces give only a mean-field contribution through 
the term C(0). Thus, the two parameters that dictate the evolution of the stochastically 
forced HMF model are C(0) and ci. From ([6]), we know that C(0) = 4:{k,)ss is proportional to 
the kinetic energy in the final stationary state. Moreover, Eq. (jl]) implies that ci < C(0)/2. 

If Ci = 0, the kinetic equation reduces to a linear Fokker-Planck equation with diffusion 
coefficient C(0)/2. This equation also describes the HMF model coupled to a Langevin 
thermostat, studied in [23l[2^ . As the kinetic equations are the same, the dynamics coincide 
at leading order in a. However, we know that at higher orders, detailed balance is broken 
in our case, whereas it holds for the Langevin dynamics. 




(13) 
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In the case Ci = 0, the homogeneous stationary states of the kinetic equation have 
Gaussian momentum distribution f{p). As has been studied thoroughly in the context of 
canonical equilibrium of the HMF model, these states are stable for kinetic energies greater 
than 1/4, i.e., for C(0) > 1. 

For values of C(0) and ci such that C(0) > 1 and ci ^ C'(O), we then expect the 
stationary states to be close to homogeneous states with Gaussian momentum, so that the 
numerical iterative scheme to locate stationary states of the kinetic equation is expected to 
converge for well-chosen initial conditions. We have checked the convergence for the set of 
values of ci used in the simulations reported in the paper. 

To check the theory, we have performed numerical simulations of the stochastically 
forced HMF model. In Fig. [H we show the evolution of the kinetic energy and (p^) = 
(l/N) "^f^iPi, and compare them with theoretical predictions. In the latter case, we have 
compared the long-time asymptotic value with the kinetic theory prediction for the stationary 
state, computed using the iterative scheme. In both cases, we observe a very good agreement 
between the theory and simulations. 

For a more accurate comparison, we have obtained the stationary momentum 
distribution from both A^-body simulations and the numerical iterative scheme. The 
comparison between the two is shown in Fig. [21 left panel, where we also show the Gaussian 
distribution with the same kinetic energy. The agreement between theory and simulations 
is excellent. 




Momentum p Momentum p 

Figure 2. The figure on the left shows tlie stationary momentum distribution f{p) for 
a = 0.01, C(0) = 1.5, and ci = 0.75. The data denoted by crosses are results of A^- 
body simulations of the stochastically forced HMF model with N = 10000, while the black 
broken line refers to the theoretical prediction from the kinetic theory. For comparison, the 
red continuous line shows the Gaussian distribution with the same kinetic energy (stationary 
state at C(0) = 1.5, ci = 0). The figure on the right shows the diffusion coefficient D[f]{p) 
for the stationary momentum distribution f{p) for different values of C(0) and ci. 
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5. Conclusions 



In this work, we studied the effect of external stochastic fields on Hamiltonian long-range 
interacting systems by generalizing the kinetic theory of isolated long-range systems. Our 
theoretical results are general, being applicable to any long-range inter-particle potential, 
space dimensions and boundary conditions. In this paper, we demonstrated an excellent 
agreement between the theory and numerical simulations for one representative case. 

Here, we discussed the kinetic theory in the limit Na ^ 1. The extension to general 
values of Na is straightforward: Because of the linearity of the equations of the BBGKY 
hierarchy, the finite- and stochastic effects give independent contributions. The kinetic 
equation at leading order of both stochastic and finite-size effects is 

^ = K[f]+LM. (14) 

where is the operator described in Eq. ffTOj) and (of order 1/N) is the Lenard-Balescu 
operator [TT]. For instance, in the case Na <^ 1 and in dimensions greater than one, the 
operator L^r is responsible for the relaxation to Boltzmann equilibrium after a timescale 
of order A^, whereas the smaller effect of selects the actual temperature after a longer 
timescale of order 1/a. 

We note that an equivalent approach to derive the kinetic theory is to write an evolution 
equation for the noise-averaged empirical density g, t) = (l/N) J2^=i{^iQii'^) ~Q)^iPi{^) ~ 
p)), by analogy with the Klimontovich approach for isolated systems. The noise appears in 
the resulting equation as a multiplicative term. This equation can be treated perturbatively, 
and may be shown to lead to the kinetic equation ffTOj) . 

Let us mention some open issues. For technical simplicity, we assumed a homogeneous 
state in our approach. Recently, Heyvaerts [25] has generalized the Lenard-Balescu equation 
to some non- homogeneous cases; his approach could be used to generalize the theory 
developed here to inhomogeneous states. There is no difficulty in principle, although actual 
computation could be more involved. 

An interesting follow up of this work is to study the dynamics of the kinetic 
equation ( fTOl) . both analytically and numerically. This may unveil very interesting behaviors, 
such as bistability or limit cycles. Bistability was observed in two-dimensional turbulence 
with stochastic forcing [26], in a framework which has deep connection with the one studied 
in this Letter. One of the motivations for this work was to make a first step in formulating 
a kinetic theory for the point vortex model and the Euler equations in two-dimensional 
turbulence [8]. This subject will be the topic of further investigations. 
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